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A time-reversal invariant Kitaev-type model is introduced in which spins (Dirac matrices) on the 
square lattice interact via anisotropic nearest-neighbor and next-nearest-neighbor exchange inter- 
actions. The model is exactly solved by mapping it onto a tight-binding model of free Majorana 
fermions coupled with static Z2 gauge fields. The Majorana fermion model can be viewed as a model 
of time-reversal invariant superconductor and is classified as a member of symmetry class Dili in 
the Altland-Zirnbauer classification. The ground-state phase diagram has two topologically distinct 
gapped phases which are distinguished by a Z2 topological invariant. The topologically nontriv- 
ial phase supports both a Kramers' pair of gapless Majorana edge modes at the boundary and a 
Kramers' pair of zero-energy Majorana states bound to a 0-flux vortex in the vr-fiux background. 
Power-law decaying correlation functions of spins along the edge are obtained by taking the gapless 
Majorana edge modes into account. The model is also defined on the one-dimension ladder, in which 
case again the ground-state phase diagram has Z2 trivial and non-trivial phases. 

PACS numbers: 75.10.Jm, 73.43.-f, 75.10.Kt 



I. INTRODUCTION 

Topological phases are a gapped state of matter which 
does not fall into a conventional characterization of con- 
densed matter systems in terms of symmetry breaking. 
The prime and classic example is the fractional quantum 
Hall effect which is realized in two-dimensional electron 
gas under strong magnetic field. Topological phases in 
the fractional quantum Hall effect are characterized, e.g., 
by the presence of (chiral) edge states, by a set of frac- 
tionally charged quasiparticles which obey fractional or 
non-Abelian statistics, and also by the topological ground 
state degeneracy when a systemjs put on a spatial man- 
ifold with non-trivial topology.El A fractional quantum 
Hall state cannot be adiabatically deformed into a trivial 
state of matter such as an ordinary band insulator. 

While it is necessary to break time-reversal symme- 
try (TRS) to realize the fractional quantum Hall effect, 
a topological phase can exist without breaking TRS, as 
seen in several examples of gapped quantum spin liquid 
states (e.g., Z2 spin liquid states). Furthermore, a phase 
which is not topological, in the sense that it can be adia- 
batically connected to a trivial phase (vacuum), can still 
be topologically distinct from the vacuum once we impose 
some discrete symmetries, such as TRS; such phases can 
be called symmetry protected topological phase.Brcl 

Symmetry protected topological phases are recently 
realized in the discovery of non-interacting topological 
band insulators, such as the quantum spin HaJL, effect 
and the three-dimensional topological insulator;&tl If we 
enforce TRS, these band insulators cannot be adiabat- 
ically connected to a trivial band insulator, as seen 
from the presence of edge or surface states. Phases of 
non-interacting fermion systems (including Bogoliubov- 
de Genne quasiparticles in the presence of meanfield BCS 



pairing gap) have been fully classified in terms of pres- 
ence or absence of discrete syHMnetries of various kind for 
arbitrary spatial dimensions .LTeI 

Studies on realizations of strongly interacting counter- 
parts of these time-reversal symmetric topological band 
insulators, i.e., "the fractional topological insulator," are 
still in their early stage.EHl 

The notion of symmetry protected topological phases 
is not limited to electron systems with TRS, but agplies 
to bosonic systems including quantum spin systems.H The 
Haldane phase in integer spin chains has been known as 
an example of a gapped spin liquid phase in one spatial 
dimension with a localized end state which carries half- 
integer spin. It is recently uncovered that the HaLdane 
phase has a symmetry protected topological order .Bu 

The list of experimentally established realizations of 
strongly interacting topological phases is still limited. 
However, a number of exactly solvable models have been 
proposed, helping us to deepen our understanding of the 
topological orders in many-body systems. Examples are 
the AfBeck-Kennedy-LiebjTasaki (AKLT) modeljO the 
quantum dimer models, 13 the toric code model,ll3 and 
the string-net models,Ej etc. In Ref. |l5|, Kitaev intro- 
duced an exactly solvable quantum spin model on the 
two-dimensional honeycomb lattice. A central feature of 
the honeycomb lattice Kitaev model, among others, is 
that it realizes, in the absence of TRS, a gapped phase 
with a chiral Majorana edge state, and non-Abelian any- 
onic excitations in the bulk. VariantS|_n£_the Kitaev 
model, such as ^jJ(2) invariant models,t3'EZl have been 
studied recently.li3 

In this paper, we consider an extension of the Ki- 
taev model on the square lattice that respects a TRS 
of some sort. Following similar extensions of the 
taev model on the three-dimensional diamond latticd^^ 
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and on the two-dimensional square lattice,Ell we consider 
two spin- 1/2 degrees of freedom on each site that com- 
pose 4x4 Dirac matrices (7 matrices). Similarly to the 
original Kitaev model, we consider interactions among 
spins which are anisotropic in space and are designed in 
such a way that the model is solvable through the (Majo- 
rana) fermion representation of spins. Written in terms 
of the fermions, our model belongs to the symmetry 
class DIILin the Altland-Zirnbauer classification of free 
fermions. LrE3 Symmetry class Dili is a class of fermions 
which are subjected to TRS, and also to particle-hole 
symmetry [i.e., a Majorana (or real) condition]. This 
should be contrasted with the original Kitaev model, 
which when rewritten in terms of Majorana fermions, be- 
longs to symmetry class D, which is a class of Majorana 
(real) fermions without TRS. One of our main findings is 
a topological phase which is characterized by the Z2 topo- 
logical invariant of class Dili in the bulk, and supports 
gapless non-chiral Majorana fermion edge modes which 
form a Kramers pair: the Bloch wavefunctions of the 
"emergent" Majorana fermions in this phase are in the 
same topological class as those of fermionic quasiparti- 
cles in the topological superconductor in symmetry class 
DHL This phase is a time-reversal symmetric analog of 
the non-Abelian phase of the honeycomb lattice Kitaev 
model, and in fact, the model can be viewed as a "dou- 
bled" version of the original Kitaev model; just like the 
quantum spin Hall system with non-trivial Z2 topologi- 
cal invariant can be constructed from two copies of the 
integer quantum Hall systems with opposite chiralities. 
Prom this point of view, our model is somewhat analo- 
gous to time-reversal invariant "doubled" anyon models 
discussed in Ref. Esl. 




FIG. 1. Square lattice and link vectors with fi — 0,1,2,3 
emanating from a site on the A-sublattice (open circle) to a 
neighboring site on the B-sublattice (filled circle). The dashed 
lines indicate a unit cell. 

II. MODEL 

In this section, we introduce an extension of the Ki- 
taev model that respects time-reversal symmetry. The 
Hamiltonian is written in terms of Dirac matrices de- 
fined on each site of the two-dimensional square lattice. 
We first consider the Hamiltonian with nearest-neighbor 
couplings only and show that its ground-state phase di- 
agram has a gapped phase and a gapless phase. O We 
then add next-nearest-neighbor couplings to the Hamil- 
tonian, which open a gap in the gapless phase. This 
gapped phase can be viewed as a topological supercon- 
ducting phase when the Hamiltonian is transformed to 
a free Majorana tight-binding Hamiltonian. The time 
reversal symmetry is preserved in both of the gapped 
phases. 



This paper is organized as follows. In Sec. ||, the 
Hamiltonian with nearest-neighbor and next-nearest- 
neighbor interactions is presented in terms of Dirac ma- 
trices and transformed to free Majorana Hamiltonian 
that respects TRS. The symmetry class in the Altland- 
Zirnbauer classification is specified, and the phase dia- 
gram of the ground states is obtained. In Sec. HI, we 



show by numerical calculation and a topological argu- 
ment that the helical Majorana edge modes appear in 
the phase with a nontrivial Z2 topological invariant. In 
Sec. 1^ some spin correlation functions are calculated 
along the edge. The existence of the gapless Majorana 
edge modes determines the power-law decay of the cor- 
relation functions. In Sec. M, we confirm that an iso- 
lated vortex excitation of the Z2 gauge field hosts a time- 
reversal pair of zero-energy Majorana bound states. In 
Sec. VI, we study the model on one-dimensional lattice. 
Two distinct phases are found which are characterized 
by the Z2 topological invariant. In the Appendix we give 
an alternative representation of Dirac matrices in terms 
of Jordan- Wigner fermions which keeps the same four- 
dimensional Hilbert space at each site. 



A. Hamiltonian with nearest-neighbor interactions 

only 

There are a class of exactly solvable quantum spin 
models in which Ising-type nearest-neighbor exchange in- 
teractions have different easy-axis directions for each link 
on the lattice. Ja the original Kitaev model on the hon- 
eycomb latticejl3 three components of the Pauli matrices 
are assigned to the three links emanating from a site of 
the honeycomb lattice. Similarly, to define an exactly 
solvable spin model on the square lattice, we can take 
Dirac matrices and assign four components of the Dirac 
matrices to four distinct types of links that emanate from 
each site, as in the Kitaev-type model on the diamond 
lattice.y 

The sites on the square lattice are divided into A- and 
B-sublattices. Pour links from a site on the A-sublattice 
are labeled, respectively, by /j, = 0,1,2,3 counterclock- 
wise from the positive a;-direction (Pig. Taking the 
lattice constant ag = 1, four types of link vectors are 
written in the two-dimensional coordinate as 



Co = 



, ei = 



62 = 



63 = 







(1) 
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where the direction of each vector is chosen from a site of 
the A-sublattice to a neighboring site of the B-sublattice. 
In the following the links labeled by 0,1,2,3) are 
referred to as "/i-links" . 

For each site on the square lattice, we consider a four- 
dimensional bosonic Hilbert space. The four-dimensional 
Hilbert space can be considered as that of a spin-3/2 
operator ,EJ or the direct product of spin-1/2 degrees of 
freedom and two orbital degrees of freedom. To describe 
the local bosonic Hilbert space, we define a set of Dirac 
matrices-p in terms of the 7 matrices in the standard 
manner :lj 



a" =7", a°=7V (a =1,2,3). 



(2) 



The Dirac matrices a satisfy the anticommutation re- 
lations {a^,a'^} — 26^'^, while the 7 matrices satisfy 
{7/^,7'^} = 23^", where g^"' = diag(l, -1, -1, -1). With 
the fifth component of the 7 matrices, 7^ = i7°7^7^7'^, 
we define another set of Dirac matrices C, 



C°=7', C=7V (a =1,2,3). 



The model is also invariant under a kind of time- 
reversal symmetry operation which is designed to become 
a time-reversal symmetry operation for half-integer spin 
fermions in the Majorana representation discussed later. 
Let us first consider a time-reversal operation T defined 
by 

T = {ia"^) (g) {iT'^)IC, 



(9) 



with complex conjugation operator )C and a = 1,2,3. 
Note that = +1. Under T, a and ( are transformed 



(3) 



as 

Tij^-f"T-^ = -^7^7°, (10) 

where covariant and contravariant vectors are defined as 
= (a*',a'^) and = (a",— a°). As we have noted, 
while the cr-part of our Hamiltonian is fully anisotropic in 
a space, the r-part of the Hamiltonian is invariant under 
a rotation around axis. In particular, it is invariant 
under a rotation i? by 7r/2 around axis. 



The C matrices also satisfy the anticommutation relations 
{C'',^^} = 25'^'^. We represent the 7 matrices as the 
direct product of two Pauli matrices cr* and (the Dirac 
representation) , 

7"=cr0(g)r3, = ia'' (g) (a = 1,2, 3), (4) 

where (T° and r*' are 2x2 unit matrices. The two sets of 
Dirac matrices are then written as 

a^' : a° = a° g) , a" = a" ® (a = 1,2,3), (5) 
: (0 = (T°®r\ C = -<^'"®T^ (a = 1,2,3). (6) 

We introduce the nearest-neighbor spin Hamiltonian, 



7^0 



3 

E 

/X— /j,-links 



(7) 



where is the coupling constant on /i-links. The sub- 
scripts j and k refer to nearest-neighbor sites on the A- 
and B-sublattice, respectively, which are connected by a 
/i-link. That is, the position vectors of the sites j and fc, 
Tj and Tfc, are related by = rj + e^. Without loss of 
generality, we can assume > Q- In terms of the two 
Pauli matrices and t^, the model can be written as 

^0 = - E -^^^ E + (8) 

/^— /i-links 

While the part of the exchange term involving the a- 
matrices is anisotropic, the part involving the r-matrices 
is isotropic and XY like; the model has a U(l) symmetry 
rotating the r-matrices around the axis. This U(l) 
symmetry, however, will be lost when we later perturb 
the nearest neighbor model (0). 



R 



R = 



V2 



(11) 



T I \ — T 

Under R, a and C arc transformed as 

Rai'R-^ = -C RCR-^ = +a'', 



Rij^-f°R-^ = +^7^7*^. 



5„,0 



(12) 



By combining T with R we can define yet another antiu- 
nitary operation, T' = i?T, 

r ^ RT ^ ^{iT^ ~ T")ia^lC, 
V2 



(13) 



Below, with a slight abuse of language, we will call this 
operation T' time-reversal operation. When applied to a 
and 



(14) 



i.e., time-reversal operation T' exchanges a and C, and 
covariant and contravariant vectors. Notice that 



r'2 = iT^ r'4 = -1. 



(15) 



We will impose the time-reversal symmetry T' through- 
out the paper. 

The Hamiltonian (^ has the integrals of motion de- 
fined for each plaquette p, 

W,^ U a^al^ n CrC' (16) 

where (j, k) are the four links on the boundary of a pla- 
quette p, and the sites j and k are on the A- and B- 
sublattices, respectively. 
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B. Mapping to Majorana fermion model 



The honeycomb lattice Kitaev model can be mapped 
to a Majorana fermion problem in the presence of a Z2 
gauge field by representing the Paiili, matrices in terms 
of four Majorana fermions per site. 113 Similarly, we can 
represent the two sets of Dirac matrices Mi^i-^^ with 
six Majorana fermions (p = 0, • • • , 5):Ej'E 



iA'^A^ 



(17) 



where we have not written the site indices explicitly. The 
Majorana fermions satisfy (A'')''' = A^ and {A'', A^ } = 
2gpp rpj^g bosonic Hamiltonian is then mapped to, 
by using the relation (17), a Majorana Hamiltonian 



where 



"Ho 



,1^ - 



(18) 



/i— /^-links 



lA^ A^^ are defined on the /i-link connecting 



two neighboring sites j and k which belong to the A- 
and B-sublattices, respectively. We will use simplified 
notation Ujk for since is uniquely determined by 
the neighboring sites j and k. The identity [uj^f' = 1 
implies that the eigenvalue of Ujk takes ±1. The Ujk 
defined on each link of the square lattice are Z2 gauge 
fields. 

Since Ujk commute with each other and also with the 
Hamiltonian (^8|) , all m^j. and the Hamiltonian can be di- 
agonalized simultaneously. Hence the total Hilbert space 
C for Majorana fermions is decomposed into subspaces 
Cs^ujk} which are specified by the configurations of the 
eigenvalues of Z2 gauge fields Ujk on every link. 



C 



3£ 



(19) 



Within each subspace, the Hamiltonian is regarded as a 
free Majorana fermion Hamiltonian, where Ujj. are re- 
placed by their eigenvalue ±1. r— ■ 

According to Lieb's theorem,c3 the energy of the free 
Majorana Hamiltonian ( p^ ) is minimized when Z2 gauge 
fields Ujk are such that each plaquette has a 7r-flux, 



n 



-1. 



(20) 



The left-hand side of Eq. (|2C|), which we denote Wp, is 
the Majorana fermion representation of the plaquette op- 
erator Wp in Eq. ( p^ and is Z2 gauge invariant. The con- 
dition (EQ) is satisfied, for example, by setting Ujk = — 1 
on the 0-links and Ujk = +1 on the other links. How- 
ever, there is redundancy in the choice of Z2 gauge-field 
configuration for a given flux configuration. 

The time-reversal operation for Dirac matrices [Eq. 
(|r^)] is translated into that for Majorana fermions as 



X 



A^ 



(21b) 



In order to keep the Z2 gauge operators invariant under 
time-reversal transformation, we employ the two types 
of time-reversal rules to M ajora na fermions on each sub- 
latti ce se parately, i.e., Eq. ( 21a ) for the A-sublattice and 
Eq. (Elbl) for the B-sublattice. 



C. Projection 

The Majorana fermion representation preserves 
the commutation and anticommutation relations of the 
Dirac matrices a and ^. However, on each site, the origi- 
nal four-dimensional Hilbert space is doubled in the Ma- 
jorana fermion representation which employs six flavors 
of Majorana fermions (or, equivalently, three complex 
fermions), as in the original Kitaev model.ll3 This redun- 
dancy can be removed by imposing the condition at every 
site I on the square lattice, 



(22) 



p=0 



The operator Di is the Majorana fermion representation 
of ililllilfli that is a unit matrix by deflnition of the 
7 matrices. The condition (p^ ) is implemented by the 
projection operator 



(23) 



acting on the states of the Majorana Hamiltonian. (In 
the Appendix an alternative representation of Dirac ma- 
trices is given in terms of Jordan- Wigner fermions which 
are free from the redundancy.) 

Now we show that the projection operator (|2^) elimi- 
nates the arbitrariness of the choice of the Z2 gauge fleld 
for a given flux configuration jWp}. Let \^]{ujk}) be 
an eigenstate of Hamiltonian (18) with a Z2 gauge-field 
configuration {ujfc}. It follows from the relation 



[Ho, A] = [M^p,A] =0, 



(24) 



that Di\^;{ujk}) is also an eigenstate of 'Hq with the 
same flux configuration {W^p}, but with a different Z2 
gauge-field configuration where the Z2 gauge fields on 
the four links around the site I are multiplied by — 1. 
This can be seen from the relations 



(21a) 



{ujk,Dj} ^ {ujk,Dk} = 0, 
[ujk,Di]=0 {l^j,k). 



(25) 
(26) 
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Furthermore, we can consider states generated by acting 
Di on multiple sites, 

nA|*;Ka), (27) 

les 

where S' is a set of sites from the square lattice. One 
might think that the total number of such states is 2^*°' , 
where A^tot is the total number of the lattice sites, since 
{Di)^ — 1. However, under the periodic boundary condi- 
tion, the number of Z2 gauge- field configurations {ujk} 
generated in this way turns out to be 2^'°*^^, since the 
product of Di on the all sites, 

X{Di ^\{u,uX{i\t\l (28) 

; (jk) I 

does not change the Z2 gauge-field configurations. More- 
over, eigenstates of the free Majorana fermion Hamil- 
tonian are invariant under the action of (^8|) up to an 
overall sign, as creation/annihilation operators of single- 
particle states anticommute with (p8|). Obviously, the 
states generated by acting Di from distinct sets S and S' 
are orthogonal, 

(*;Kfc}in^' n -0, (29) 

leS I'es' 

unless S* = S" or 5" is the complementary set of S, since 
the eigenvalues of the Z2 gauge-field operators are dif- 
ferent between two states. Hence, the states of the form 
( p7| ) form 2^'°' ^^-dimensional orthonormal basis states. 
Meanwhile, the number of flux configurations is 2^'°'~^, 
since the total fiux must be unity (J^^ Wp — 1). Consid- 
ering the fact that there are two additional, independent 
integrals of motion defined on two closed loops Cx,Cy 
going around in the x- and y-directions, 

W^^ n "X' (30a) 

(j,fc)ec, 

Wy= l[ a^^al (30b) 

we find that the number of Z2 gauge-field configurations 
for a given local fiux configuration ({Wp}) and global 
flux conflgurations {{Wx,Wy}\]s 22^*°V(2^'°*-i22) = 
Therefore the states ( p7[ ) exhaust the eigenstates 
for all Z2 gauge-field configurations with the same flux 
configuration. Finally, projecting the state ( p7| ) yields 

pn^n*;K4) (3i) 

les 

since {1 -\- Dj)Dj = 1 + Dj. Equation (|3l| ) implies that 
the projected state is independent of Z2 gauge choice. 
Whatever Z2 gauge configuration is taken for a given 
flux conflguration, the same set of states are obtained 



after the projection; any redundant state of the free Ma- 
jorana fermion Hamiltonian disappears after the projec- 
tion. Furthermore, we can conclude that matrix ele- 
ments (for the projected states) of gauge-invariant ob- 
servables can be calculated by using eigenstates of Ma- 
jorana fermions with any particular Z2 gauge configura- 
tion. 



D. Phase diagram of the nearest-neighbor spin 
Hamiltonian 

Let us set Ujk = — 1 on every 0-link and = +1 on 
the other links of the square lattice, to satisfy the 7r-flux 
condition, Eq. (pO|). This Z2 gauge-field configuration, 
which we denote by a four-vector m'' = (m°, u^, u^, u^) = 
(—1,1,1,1), preserves lattice translation symmetry with 
the unit cell shown in Fig. ||. We introduce Fourier 
transformation of Majorana fermion operators on the A- 
sublattice, 

= -L= ^ e-'-^-^ (s = 4, 5), (32a) 

and of those on the B-sublattice, 

6^ = -i=^e-''-^A^ (s = 4,5), (32b) 

where N is the number of unit cells, and rj in both of 
Eqs. ( p2| ) is the position vector of th e site j on the A- 
sublattice. That is, Vj in Eq. ( ^2b| ) is related to the 
position vector of the site k on the B-sublattice by r,- = 
rk + ei. The inverse Fourier transform of Eqs. ( |3^ ) is 
given by 

q 

KeB ^K = \[^T.^"'-'''"'^'''K (33b) 
q 

where the wave vector q is in the first Brillouin zone. 




I + kyl ^ The fermion operators defined in Eqs. 
) satisfy the following relations: 



a% = (a^)\ (34a) 

(34b) 

One can thus regard and alg (6^ and folg) as annihi- 
lation and creation operators of fermions (or vice versa) . 
Hamiltonian ( p^ ) is written in the momentum space as 

Ho-5][i'i>(q) {a%h\ + at^bl) 
q 

-*$*(q)(5%4 + 6l,a^)], (35) 
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(0,0,0,1) 



(1,0,0,0) 




(0,1,0,0) 



FIG. 2. Phase diagram in the parameter space (Jo, Ji, J2, J3). 
Shaded regions are gapped phases and the other area is gap- 
less phase. 



where 

$(q) = e'"" J^u^e*«-^f' 



(36) 



The eigenvalues of ( |35| ) are i? = ±|<i>(q)|. Each eigenstate 
is doubly degenerate, since and are decoupled in 
the Hamiltonian. The ground state is obtained by filling 
all the eigenstates with negative energy. 

The ground-state phase diagram is drawn in Fig. |^, 
where for illustration purpose we normalized the param- 
eters J = (Jo, Ji, J2, J3) such that Jo + Ji + J2 + = 1, 
Jfi > 0. The vertices of the (large) tetrahedron in 
Fig. I, J = (1,0, 0,0), (0,1, 0,0), (0,0, 1,0), (0,0, 0,1), 
correspond to the parameter sets in which one of the four 
coupling constants is much stronger than the others. On 
the edges of the tetrahedron the sum of two coupling 
constants is equal to 1. The four shaded regions (smaller 
tetrahedrons) in Fig. ^ are gapped phases in which there 
is an energy gap between positive energy bands and neg- 
ative energy bands. The region including the isotropic 
point J = (1/4,1/4,1/4,1/4) (the non-shaded part in 
Fig. H) is a gapless phase where the positive and negative 
energy bands touch at two Dirac points, around which 
Majorana fermions have linear energy dispersions. The 
gapless phase will become a gapped topological phase, 
once an energy gap is opened by some perturbations, as 
is the case in the Kitaev model. On the boundary be- 
tween a gapped phase and the gapless phase, two Dirac 
points merge to become a single point in the Brillouin 
zone. This happens when one of the four J^'s is equal to 
the sum of the other three J^'s. 



E. Hamiltonian with next-nearest-neighbor 
interaction 

We add, to the Hamiltonian Ho, perturbations of the 
form of a product of Dirac matrices from three neighbor- 



ing sites. As we will see, these perturbations will open a 
gap at the Dirac points in the gapless phase. 

Consider three neighboring sites j, k, and I of a single 
plaquette shown in Fig. |^, where the sites j and k belong 
to the same sublattice (cither A or B). We consider three- 
site interaction Hamiltonian of the form 



(37) 



where the links (j7) and {Ik) are a /i-link and a z/-link, 
respectively. In the Majorana fermion representation, the 
three-site interactions read as 



i(a^af)(«) = m^,<;A|At 

»(c;cr)(cra) = K/<iA,5At 



(38a) 
(38b) 



As the Majorana operators Af do not appear explicitly 
in the right hand side of Eqs. (38), we can regard these 
perturbations as next-nearest-neighbor hopping opera- 
tors for A^'^. We have a different type of three-site inter- 
actions in which different sets of Dirac matrices are used 
for two links: 

dm 

(39) 

This yields another type of next-nearest-neighbor hop- 
ping term. 



*(cj'cr)7f7f («r«D = *<,"fe/A,5At 



(40a) 
(40b) 



The summations in Eqs. (^7|) and ( |39D are over any three 
neighboring sites which belong to the same plaquette, 
including the three sites (jl'k), in addition to {jlk), in 
Fig. ||. For the Z2 gauge fields satisfying the tt-Aux con- 
dition, Eq. (|o|), the product of Z2 gauge fields UjiUki in 
Eqs. ( ^Sf ) has the opposite sign compared to the product 
Uji>Uki'. We assume Js:|,^ = -iiTJ;,^ and K^i,^ = -K'^vk 
so that the two paths give the same contributions. This 
leads to a vanishing next-nearest-neighbor hopping for 
0-flux plaquettes. 

To summarize, under the tt-Aux condition ([20|), we 
have the following two types of next-nearest-neighbor 
hopping Hamiltonian: 



jGA a=l,2 

■E E 

■Y, Y ^^:(A|Af+„„ 

jSA a=l,2 

-E E ^^(^^^.. 

j£B a=l,2 



AjAj+Q„) 



AjAj+a„), 



AjA|+„^) 



(41a) 



+ >^^+aJ, (41b) 



where the parameters i^^'^ are hopping matrix elements. 
The subscript j -t- denotes the site located at Tj + aa 
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0k 



FIG. 3. Four-spin interaction terms. Left: The Z2 gauge 
fields obtained from the two paths connecting next-nearest- 
neighbor sites. Right: Two directions are labeled by a = 1, 2, 
respectively. 



with ai = (1, 1) and 02 = (—1, !)• The two vectors 01^2 
correspond to the two directions of next-nearest-neighbor 
hopping labelled by a = 1 , 2 in Fig. |^. We shall see that 
the next-nearest-neighbor hopping terms of the form 

and z(A^A| + A^At) (42) 



Here we have defined d {i = x, y, z) as the Pauli matrices 
acting on the sublattice indices (a, h), and as those on 
the Majorana flavors (4,5). The matrix s° is the 2x2 
unit matrix in the Majorana flavor space. The Hamilto- 
nian in the momentum space x{q) is invariant under the 
translation by reciprocal lattice vectors, G = (±7r,7r). 
The eigenenergies are ±£q, where 



(48) 



Each energy level, for a given q has two-fold degeneracy. 

We have seen in Sec. 11 D that the gapless phase of 
Hamiltonian has two Dirac points at q — Qi where 
$(qi) = 0. Non-vanishing matrix elements &x{q) and 
Qz{q) at the Dirac points give a band gap. Hence the 
gapless phase is turned into a gapped phase by including 
the three-spin interaction or next-nearest-neighbor hop- 
ping interactions Hz and Hx- 

Incidentally, both Qx{q) and &z{q) vanish on the 
phase boundaries of the gapless and gapped phases of 
Hq. Thus, the phase boundaries do not change upon 
addition of Hz.x to Ho- 



are invariant under time-reversal transformation which is 
defined in the next section. 



Fourier transformation of Eqs. (41a) and ( fflq ) yields 

,4 l4 .4 „5 „5 I l5 1,5 ^ 



Hz =5]e,(q)(ai,c 
=5]e,(q)(a%c 



(43a) 
(43b) 



where 



e,{q) = Kl sin{qx + Qy) + Kf sin(-g^ + qy) (44) 
for i — z,x. The total Hamiltonian can be written as 



n = no + nz + nx = J2 ^U{q)%, 



(45) 



where 



if 

\4j 



(46) 



and 



V 



i$ 

-Qz -e, 

-e^ i<p 

-Ox 



= - Re [$(q)] (g) s° - Im [$(q)] , 
+ Qziq) (g Qxiq) 



(47) 



F. Symmetries 

In this subsection, we consider symmetry properties of 
the Hamiltonian H in the Majorana fermion representa- 
tion and show that it belongs to DILL symmetry class of 
the Altland-Zirnbauer classification.CJ 

To this end, we begin with transforming the Majo- 
rana Hamiltonian H into Bogoliubov-de Gennes (BdG) 
Hamiltonian. We define fermion creation and annihila- 
tion operators on site j from the two flavors of Majorana 
fermions A| and A^, 



1 



Ci 



1 



2' J -J" -3 2 
Their Fourier transforms are written as 
1 



(49) 



4, 



E 



for the A sublattice (j € A), and 

1 ■r-^ fg-iq-irk+ei) 



Bq ^ _ 

Bi 



Ck 

oiq-(rk+ei ) J 



(50a) 



(50b) 



for the B sublattice (fc G B). The Nambu field for the 
complex fermion is defined by 



(Aq\ 

Bq 

A-q 



(51) 



and is related to ipq in Eq. (| 
mation, 

i^q = C/*, 



by the unitary transfor- 



(52) 
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where the unitary matrix U is given by 



where 




(53) 



Then the Hamiltonian "H can be written in the form of 
BdG Hamihonian, 



(54) 



where 





e*(q) 










-e*iq) -i<f*iq) 



Re [$(q)]c^ 
Re [e(q)]c^ 




Im [$(q)]c'= (g> t" 
-lm[Q{q)]c' (E)ty 



with 8(q) defined by 

e(q) = e,(q) + ie,(q). 



(55) 



(56) 



(57) 



The Pauh matrices {i = x, y, z) and the 2x2 unit 
matrix t'^ act on the Nambu indices. 

We are ready to discuss symmetries of our model in 
terms of the BdG Hamiltonian x{q). Under the particle- 
hole transformation generated by T' = t^K., where JC 
is complex conjugation operator, the BdG Hamiltonian 
changes its sign, 



(58) 



We note that "P^ = +1. 

The BdG Hamiltonian is invariant under time-reversal 
transformation 



i{c'®ty)x' {~q){~iW®ty)^x{q)- 



(59) 



The time-reversal operator T = c^^ityjC obeys — —1. 

We conclude from these symmetry properties that the 
BdG Hamiltonian x belongs to symmetry class DHI; see, 
e.g.. Table 1 in Ref. ^. It is known from the classification, 
theory of topological insulators and superconductorsLrcl 
that gapped ground states of class DHI Hamiltonian in 
two spatial dimensions can be classified by a Z2 index; 
see, e.g.. Table 3 in Ref. ^. 

The product of particle-hole and time-reversal trans- 
formations, TV ^ yields 



(60) 



i.e., x{q) anticommutes with (g) i^. In the basis where 
c^(g)t^ is diag(l, 1,-1,-1), the BdG Hamiltonian is writ- 
ten in the off-diagonal form. 



XD{q) = l2iX{q)hi 



( D{q) 
Wiq) 



(61) 



hi — 



a o\ 
0001 
0010 
^0 1 oy 



(62) 



and 



D{q) = 



Q{q) t^q) 
-t^*{q) -e*{q) 



(63) 



Since I2A and commute, we find from Eqs. ( |58| ) and 
( |6l| ) that D{q) satisfies the skew relation 



D'^i^q) = -D{q). 



(64) 



Since /24C^ ® = it^ , the time-reversal operator for 

XD{q) is T=ityiC. 



The symmetry relations in Eqs. (^8|) and (59) lead to 
symmetry relations for the Majorana Hamiltonian xiq) 
through the unitary transformation. The particle-hole 
symmetry relation implies 

X'^i-q) = -X{q), (65) 

while the time-reversal symmetry gives 

® (is^)x^(-g)c" ® i-isy) = xiq)- (66) 

It follows from these relations that xiq) anticommutes 
with (8) s**. 



(c"®s^)x(g)(c"<»s^) = -x(g). 



(67) 



We note that the time-reversal operator c^(g)(?s*')/C in Eq. 
( |66| ) is consistent with that for the Majorana operators 
[Eqs. (|l|) and (|21b|) ]. 

We return to the time-reversal symmetry of the next- 
nearest-neighbor hopping terms, Eq. (^2|). When the 
time-reversal operator T = ® (is^)/C is applied to the 
A*A^ , does nothing since both sites j and k are on the 
same sublattice, while is^ interchanges s = 4 and s = 5 
with a factor of —1 (+1) for s ^ s' [s = s'). Thus, 



(68a) 
(68b) 



and the hopping terms in Eq. (^2|) are invariant under 
the time-reversal transformation. 



III. PHASES AND TOPOLOGICAL INVARIANT 

In this section we show the existence of a Kramers' pair 
of Majorana edge modes in the topological phase and de- 
fine a Z2 index that distinguishes between the topologi- 
cally nontrivial and trivial phases. 
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(a) 
Ji/J 



(b) 



I A-phase 



-o-KS-^^^e-K)- 



(4)' 



B-phase I 



{3)' 



.(2) 



A-phase 



~9 — 9 — 9 — r 



•(1) 



■hl-J 



FIG. 4. (a) Phase diagram in the Jo/ J — Ji/J plane. The 
region between the dashed Unes Ji = Jo ± 2J, including the 
isotropic point Ji = Jo = J, is the gapless phase when the 
next-nearest-neighbor hopping terms are absent ("B-phase"). 
The gapless phase is turned into a topologically nontrivial 
phase when the second-nearest-neighbor hopping is turned 
on. Both the gapped phases for Ji > Jo + 2J and for Ji < 
Jo— 2 J remain topologically trivial upon including the second- 
nearest-neighbor hopping terms ("A-phase"). The numbered 
solid circles indicate the parameter sets for which the energy 
spectra are calculated and shown in Fig. H (b) Strip geometry 
with edges along the vector eq. 



A. Energy spectrum 

We examine the energy spectrum of Majorana fermions 
by varying two coupling constants Jq and Ji while the 
others are kept fixed as J2 = J3 = J and K], = — 
K]. = K'^ = K, for simplicity. All coupling constants are 
taken to be positive. We set the Z2 gauge fields — 
(— 1, 1, 1, 1) as before. 

When K = 0, the eigenvalues of x{q) are given by 



u=0 



(69) 



When Jq — 2J < Ji < Jq -I- 2 J the positive and negative 
energy bands touch at two Dirac points ("B-phasc" in 
Fig. y). For example, in the isotropic case Jq = Ji = J 
the Dirac points are located at q = (0,±7r/2). As we 
approach the phase boundaries Ji — Jq = ±2 J, the two 
Dirac points come closer to each other and eventually 
merge at q = (0, 0) for Ji — Jq— 2 J and at q = {tt/2, n/2) 
for Ji — Jo + 2 J. On the other hand, in the gapped phase 
where Jq — 2 J > Ji oi Jo+2J < Ji, there is an energy gap 
between positive and negative energy bands ( "A-phase" ). 

The effective Hamiltonian around the Dirac points is a 
Dirac Hamiltonian with the mass terms proportional to 
the next-nearest-neighbor hopping K. Note that 8(q) 
vanishes at q — (0,0) and (7r/2,7r/2). This means that 
the band gap is closed at Ji ~ Jq±2J even in the pres- 
ence of the next-nearest-neighbor hopping. Thus the pa- 
rameter space Jq/ J — Ji/J is divided into three regions 
by the phase boundaries Ji = Jo ± 2 J (the dashed lines 
in Fig. |). 

We have numerically diagonalized the Majorana tight- 
binding model Ho+T-Lz+l-Lx for the strip geometry, shown 



in Fig. ||(b) , where the edges are parallel to the link vec- 
tors eo and 62. The energy spectra of Ho in the strip 
geometry are shown in Fig. || for J^/ J = (4, 1, 1, 1) [(la) 
and (lb)], J^/J = (1,1,1,1) [(2a) and (2b)], J^/J - 
(3,3,1,1) [(3a) and (3b)], and J^/J = (1,4,1,1) [(4a) 
and (4b)]. We have chosen two values for the next- 
nearest-neighbor hopping: K = [(la), (2a), (3a), (4a)] 
and K/J = 0.15 [(lb), (2b), (3b), (4b)]. 

Without the next-nearest-neighbor hopping terms 
{K = 0), flat bands appear exactly at zero energy in 
the region Jq — 2J < Ji. In the gapless phase (B-phase), 
the energy spectrum shows two Dirac points at time- 
reversal symmetric momenta, and the doubly degenerate 
zero-energy flat bands connect these two points through 
Qx = (i-G-j not through Qx = 7r/2); see Fig. I (3a). In 
the gapped phase (A-phase) where Jo + 2 J < Ji, the bulk 
bands are fully gapped, and the flat bands are extended 
in the whole Brillouin zone [Fig. |^ (4a)]. 

The existence of these zero-energ]iJlat bands can be 
explained by a topological argument.EZl When K = 0, 
and decouples in our model. The bulk Hamiltonian 
for A* (or A'"^) has the form 



h{q) - Riq) ■ (T, 
where R{q) is a two-dimensional vector. 



R{q) 



-ImMq)] 
-Re[<&(q)] 



with 



$(q) = e-'i^i-Joe"^^ + Jie'iy + Jae" 



(70) 



(71) 



(72) 



q is the wave vector in the first Brillouin zone, and 
cr = {a'-'.cry). In Eq. (^2]) $(q) has the phase factor 
g-fjx j^.^ cl>(qt) in Eq. (^6D], because we have chosen the 
unit cell depicted by the dashed line in Fig. ^(b) which is 
commensurate with the presence of the boundary. Note 
that h{q) has chiral symmetry, {/i(<7),cr^} = 0. When 
we fix qx, h{q)\q^ can be regarded as a one-dimensional 
Hamiltonian with wave number qy in the direction per- 
pendicular to the edge. The one-dimensional Hamilto- 
nian has zero-energy edge modes if a loop trajectory that 
R{q)\q^ draws as qy is varied encloses the origin R ^^Jl'm. 
the two-dimensional parameter space R — (i?^, i?^)E3'E§ 
The number of zero-energy edge modes is given by the 
winding number of the loop and can change only when 
the loop touches the origin, i.e., when a band gap closes. 
For a fixed qx, the loop is an ellipse described by the 
equation 



[i?- + (Jo + J2)sinq,]2 [RV + [J^- J^)cosqx]' 



{Jl 



= 1. 

(73) 



[In deriving Eq. ([73|) we have ignored the phase factor 
in <&(g), since it does not change the winding num- 



ber.] The flat band appears when the origin i? = is 
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inside the ellipse, i.e., when satisfies 

< 4(JoJi + J^JzW, + Ji J2) ■ ^ ^ 

In the phase diagram in Fig. ^(a), the right-hand side 
of Eq. (|4|) is larger than unity for Ji > Jq + 2J (A- 
phase in the upper side), while it is less than zero for 
Ji < Jq + 2 J (A-phase in the lower side), which explains 
the spectra shown in Fig. || (la) and (4a). Otherwise, 
when Jq — 2 J < Ji < Jo + 2J (B-phase), the right-hand 
side of (jrj) takes an intermediate value between and 1, 
corresponding to the Fig. || (2a) and (3a). 

When the next-nearest-neighbor terms are included 
{K ^ 0), the bulk bands are gapped in the whole re- 
gion of the A- and B-phases. Then the fiat bands are 
split from the zero energy, except at the time-reversal in- 
variant momenta Qx = 0,7r/2. Hence, the edge modes in 
the B-phase have a single zero-energy point in the first 
Brillouin zone [Fig. I (2b) and (3b)], since the flat bands 
for K = pass through = only, while those in the A- 
phase have an even number of zero-energy points [Fig. ^ 
(lb) and (4b)]. 



B. 



index 



The phase with topologically protected edge states is 
characterized by a nontrivial Z2 index calculated iru 
bulk. The Z2 index introduced by Kane and MeldS 
for time-reversal invariant band insulators in class All is 
defined through the matrix 



(75) 



where \uj{q)) is the single-particle Bloch wave function 
in the i-th filled bands. The Z2 invariant 1^ is then given 

by 



n ^det[ui(q)] 
Pf\w(q)] ' 



(76) 



where TRIM is the time-reversal invariant momenta 
in the Brillouin zone, (0,0), (7r/2,7r/2), (0,7r), and 
(— 7r/2, 7r/2). The sign of the square root in the numera- 
tor is chosen to be continuous along the path connecting 
the four time-reversal invariant momenta. The topologi- 
cal phase of our model in class Dili can also be charac- 
terized by the Z2 index defined by Eq. (|76|). 

For models in symmetry class Dili, the Z2 index be- 
comes apparent when the Hamiltonian is expressed in 
the off-diagonal form by utilizing the chiral symmetry. 
Let us introduce the operator Q that has the eigenvalue 
-1-1 (—1) for the states in the empty (filled) band of 
the BdG Hamiltonian x- When x{q) is diagonalized as 
X = l^diag(ei, 62, • • • , — ei, —£2, • • • )V~^ with a unitary 
matrix V, the operator Q is given by 




FIG. 5. Energy spectra of the one-dimensional strip [Fig. 
^b)] as a function of the momentum along the edge. The 
energy spectra are calculated for the following parameter sets 
[see also Fig. §(a)]: J^/J = (4,1,1,1) in (la) and (lb); 
J^/J = (1, 1, 1, 1) in (2a) and (2b); JJJ = (3, 3, 1, 1) in (3a) 
and (3b); J^/J = (1,4, 1, 1) in (4a) and (4b). The Z2 gauge 
are fixed as u'^ = (—1, 1, 1, 1). The second-nearest-neighbor 
hopping = (a = 1,2 and i = z,x) in the left figures 
[(la), (2a), (3a), (4a)], K°'/J = 0.15 in the right figures [(lb), 
(2b), (3b), (4b)]. 



In the basis where x{q) takes the off-diagonal form of Eq. 
(|6l|), the operator Q also takes the form 



(78) 



Q = Fdiag(l, 1, 



-1,-1, 



)v- 



(77) 



The off-diagonal component q{q) satisfies the relations 
q^{—q) = —q{Q) [so does D{q)] and q^q = qq'^ = I (from 
— /), where / is a unit matrix. In this basis the 
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operator Q is related to the BdG Hamiltonian x by 

1 



Q{q) = —XD{q) 



(79) 



with £q defined in Eq. (^ . 

The eigenvectors of Q in Eq. (p8) are given by 



1 



(a = 1,2), (80) 



where ± indicates the eigenvalue ±1 (i.e., empty and 
filled bands), and Ua are unit vectors. 



0/ ' 



(81) 



Since the eigenspace of the operator Q of the eigenvalue 
— 1 is the same as the Hilbert space that spanned by 
the filled bands of the BdG Hamiltonian, the Z2 index 
calculated with the vectors in Eq. (|8^) is equal to that 
calculated for the original BdG Hamiltonian. 

Applying the time-reversal operator T = it^JC to the 
eigenvector of the filled states given in Eq. (|o|) yields 



(82) 



The matrix w is then obtained as 



Wabiq) = {ua-i-q)\Tub-{q)) = -qbaiq)- (83) 
It then follows from Eqs. (|l]), (jes]), and (^ that 
1 /e(q) -i<f*{qT 



w[q) 



UHq) -Q*{q) 



(84) 



Note that det[w(<7)] = — 1 in the whole Brillouin zone, 
and that w{q) becomes a purely imaginary antisymmetric 
matrix at the TRIM. Hence the Z2 index in Eq. (^^ is 
reduced to 



q:TRIM 

At the TRIM we have 

$(0,0) = -Jo + ./i + ./2 + J3, 

$(0,7r) = J0 + ./1-./2 + J3, 
$(V2,7r/2) = Jo- Ji + J2 + J3, 
$(-7r/2, V2) = - Jo - Ji - J2 + Js- 



(85) 



(86a) 
(86b) 
(86c) 
(86d) 



The isotropic point, J^ = J, has v = —1 and is thus 
a topologically nontrivial state. In the limits where one 
of J^ is much larger than the other three, v = +1 and 
the ground state is topologically trivial. These results 
are i n agre ement with the numerical results presented in 
Sec. HI A , The presence or absence of helical Majorana 
edge states is dictated by the Z2 index v. 

At the phase boundaries between the topologically 
nontrivial phase and trivial phases, ^q) vanishes at least 
at one of the TRIM. Using Eqs. (|8^), we arrive at the 
phase diagram shown in Fig. ^, in which the shaded re- 
gions are topologically trivial phases and the rest is a 
topological phase (except on the phase boundaries). 



IV. EDGE STATES AND SPIN CORRELATION 
FUNCTION 

In this section we examine spin correlations of the 
ground state of the model, especially spin correlation 
functions along the edge of the two-dimensional system. 
The edge states appear in time-reversal pairs and form 
the helical Majorana edge mod(^ as in the time reversal 
helical p-wave superconductors.Q Spin correlation func- 
tions are calculated for these helical edge states. 

Before proceeding to the calculation of correlation 
functions, we examine which operators have non- 
vanishing expectation values in the ground state. Op- 
erators can vanish due to two reasons: the projection 
operator and the constants of motion. The former one 
restricts non- vanishing operators to be the product of an 
even number of the Majorana operators on each site, be- 
cause 



1 



1 



D 



0, 



(87) 



where p = 0, . . . , 5. The latter one implies that non- 
vanishing operators can contain the Majorana operators 
with /i = 0, 1, 2, 3 only in the form that does not flip 
the Z2 flux {Wp}, since such Majorana operators alter 
the Z2 gauge configuration: 



A'; = -A'; 



3 "jfe 



(a* = 0,1,2,3), 



-l-jkXj — X'^Ujk, UjkX^ — X^Ujk, 



(89) 



and the projection operator does not flip the flux. These 
conditions can be restated as follows: non-vanishing op- 
erators for the ground state |GS) of the 7 matrix Hamil- 
tonian are 

1. Z2 gauge operators on closed strings of links, 

2. Z2 gauge operators on open strings of links, each 
string having either A^ or A^ at the both ends, 

3 . ^X'j X'j , 

4. D,, 

and products thereof. The product of the Z2 gauge field 
operators on closed strings is rewritten as the product of 
the plaquette operators Wp and gives extra minus sign, 
since the plaquette operators are the integrals of motion 
for both the Majorana Hamiltonian and the 7 matrix 
Hamiltonian. The correlation of operators which do not 
satisfy the above conditions does not extend beyond the 
nearest neighbor. 

Here we consider the two-point correlation function 
of a local spin operator, which is a single-site oper- 
ator, i.e., an operator that does not have a "string" 
composed of a product of operators. The only non- 
trivial single-site Hermitian operator that satisfies the 
above conditions is iXjXj. The term iXjXj is equal to 
^''^ — (cr° (g) r^)j in the 7 matrix and Pauli 
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matrix representation, and is also written as 2cjCj — 1 where < is an integer, 
in terms of the complex fermions (E^). In the bulk, the 
two-point correlation function of *A|A^ is short-ranged 



since A| and are free Majorana fermion operators and 
the bulk is gapped. On the other hand, the correlation 
function of iX'^X^ along the edge is expected to decay 
algebraically in the topological phase. 

Let us consider the semi-infinite system with the edge 
along the x axis, and let Qx be the momentum along the 
edge, which is conserved. The fermion operators can be 
expanded as 



t/2 



Jo 



(90) 



where r = {rx,ry) labels unit cells which contain two 
sites (the semi- infinite system is defined for ry < 0), 
fq^,i{ry) is the i-th exact single particle wavefunction 
with momentum and energy Eq^^i, and f*^ i{i~y) is 

the particle-hole conjugate with —Eq^^i\ ^^^ ^ (^^^ ■) is 
the fermion annihilation (creation) operator associated 
with these levels. In computing the correlation function 
on the edge, the dominant contributions come only from 
the modes localized at the edge. Namely, the fermion 
operators near the edge are approximated by, at low en- 
ergies, 



hi 



r 



Tr/2 



dqx [e"^-''-g+.q^{ry)j+,q^+h.c.] 



Tr/2 



+ / dq, [e*9-'--g_,,^(rj,)7_,,^ + h.c] , 
Jo 

(91) 

where g±^q^ (ry) is the single-particle wavefunction of the 
left (right)-moving edge mode with momentum q^, and 
l±,q^ is the corresponding fermion annihilation operator. 
The edge contribution to the Hamiltonian is given by 



U 



edge 



tt/2 



dqxE{qx) {-fX^q^-^+^q^ - 7I ,j^7_,,^) 



(92) 



The energy dispersion E{qx) for the edge mode is linear 
around a TRIM, g^^ = 0, as shown by the numerics in 
Fig. |. 

Ki qx =0, the edge states are doubly degenerate at 
E^Q. For Kl = Kl = and = Kl = the 

zero-energy eigen wavefunctions can be explicitly written 
as 



(93) 



^0 



^0 



{Jiu^ - J3U^)/2±A 
-Kx 
V 




(94) 



(95) 



with A = ^Kl +KI + (Jiui - J3u3)2/4, and Ai,2 are 
solutions of 



0. 



(96) 



When I All > 1 and IA2I > 1, the wavefunctions (|93|) are 
normalizable and localized near the edges. Such solutions 
of A exist when 



±A 



< 



and 



(97) 



(98) 



where all Jn are assumed to be positive. The former 
condition (B^ determines which signs to be taken. The 
latter condition (|8|) coincides with the region in Fig. ^ 
where there are zero-energy edge states at q^; = 0. In 
lowest order in q^, the two-fold degeneracy of the zero 
modes is lifted; the energy dispersion near q^ = Q is, E = 
±vq, with the velocity 



(99) 



In lowest order in q^, the eigen wavefunctions g±\^ are 
linear combinations of 5^^=o ^'^^ 9q^=o- 

From the edge theory with a linear dispersion at low 
energies, one can immediately see the (equal-time) two- 
point correlation function of the Majorana fermion op- 
erators decay along the edge as (A''(r3;)A'' (r^)) ~ {r^ — 
r^)~^. The two-point correlation function of the opera- 



tor 17^7° 

the Wick's theorem, as 



cr*^ (g) r'^ = iA^A^ can be represented, using 



C(ry,r'y) 



(100) 



where C{ry, is a function of ry and r'y which is deter- 
mined by the wavefunction of the edge modes g±^q^ and 
decays exponentially into the bulk. 
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FIG. 6. (a) The number of states in the topological phase 
for 20 X 42 system with a pair of vortices separated by 20 
lattice spacings. Periodic boundary conditions are imposed 
in the x and y directions. The parameters in the Hamiltonian 
are J*" = 1 and Ki]^ = 0.3. (b) The energy difference AE 
between the positive and negative energy eigenvalues of bound 
states, as a function of the distance between the vortices. 



V. VORTEX BOUND STATES 

In this section we discuss vortex bound states in the 
topological phase. An isolated vortex in a topological su- 
perconductor can accommodajbe a topologically protected 
zero-energy Majorana state£d The time-reversal symme- 
try of our model implies that there are two such Majorana 
zero-energy states which form a Kramers' doublet. 

In our model a vortex corresponds to a plaquette with 
a 0-flux in the 7r-flux background. Such 0-flux excitations 
always appear in pair since the total flux is a conserved 
quantity modulo 27r. As we noted above, each vortex 
should have two Majorana bound states. To confirm the 
number of bound Majorana states, we have numerically 
diagonalized Hamiltonian in the topological phase (the 
B phase in Fig. ^ for the system size of 20 x 42 sites, in 
which two 0-flux plaquettes are placed along the y direc- 
tion. We have imposed periodic boundary conditions in 
the X and y directions and set the parameters as = 1 
and K^'l = 0.3. 

Figure ^(a) shows the number of eigenstates when two 
vortices are separated by 20 lattice spacings. We find 
four nearly-zero-energy states inside the bulk gap, i.e., 
two states per vortex. The energy eigenvalues of these 
midgap states are zLAE/2, each energy eigenvalue being 
two-fold degenerate. Figure |(b) shows AE as a function 
of the distance r between the two vortices. The depen- 
dence of AE on r is symmetric about r ^ 21, because 
of the periodic boundary conditions imposed. The clear 
exponential dependence on r (r < 21) confirms that the 
energy difference AE is due to a small overlap of expo- 
nential tails of wave functions bound to the two vortices. 



VI. EXTENDED KITAEV MODEL ON THE 
ONE-DIMENSIONAL LATTICE 

In this section we study the extended Kitaev model 
on the cylinder geometry, i.e., the ladder with two sets 




FIG. 7. One-dimensional system on a two-leg cylindrical lat- 
tice. 



of rungs. The sites on the ladder are divided into A 
and B sublattices, shown as open and filled circles in 
Fig. ^ respectively. The /i-links are defined as in the 
two-dimensional case. 

We consider the nearest-neighbor interaction Hamilto- 
nian 

Hid - - E E K< + C;Cn- (101) 

fi—O /^-links 

We use the Majorana fermion representation of the Dirac 
matrices (p^), and combine and to make complex 
fcrmions as in Eq. (^9|), and take the Fourier transform 
of the Majorana operators 



jeA 



(102a) 
(102b) 



where L is the length of the cylinder, yj and yk are the po- 
sitions of the sites in the vertical coordinate. The Hamil- 
tonian in the momentum space is then given by 



Hid = * E E ^U^Pt + A,'A|) 

l-i—O /1-links 
= E*J^1D*9' 



(103) 



where the spinor ^'g is 



(104) 



and 



XlD = 



i$(q) 
-i<i>*{q) 


-i$*(Q) 






i$(g) 




(105) 
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(1,0,0) 




$(7r) = Jo + Ji + J2 - h- 



(109b) 



(0,1,0) 



(0,0,1) 



FIG. 8. Phase diagram of the one-dimensional cylindrical lat- 
tice model in the parameter space (Jo-I- J2, Ji, Js)- The phase 
in which one of the horizontal bonds Jo and J2 is greater than 
the other three bonds is a topologically trivial phase. The 
shaded regions (tetrahedra) are a topological phase, where 
each end of the ladder has two-fold degenerate zero-energy 
Majorana states. 



with 



$(<?) = Jou" + Jl 



(106) 



The eigenenergies of xiD are E — ±|$(g)| and each en- 
ergy level is doubly degenerate. From the Lieb's theorem, 
the ground state is obtained for the Z2 gauge field config- 
urations with TT-flux per each square, i.e., when the con- 
ditionsgn[(Jiui)(J3u3)(JouO + j2-u2)(Jo-uO-HJ2u2)] = -1 
is satisfied. Without loss of generality, we will work with 
the Z2 gauge = (1,-1,1,1). The ground-state energy 
is then a function of three parameters Jo + J2, Ji, and 
J3. Even without the next- nearest neighbor interaction 
terms, the ground state is gapped except at the phase 
boundaries. 



Jl - J3 = ±(Jo + J2) 



(107) 



The phase diagram is depicted in Fig. ||. 

We diagonalized numerically the Hamiltonian for a fi- 
nite length system with open boundary condition in the 
leg direction. No midgap states are found if Ji — J3 < 
Jo -t- J2 and Jl — J3 > — (Jo + J2) (non-shaded region 
in Fig. while midgap states bound to each end are 
found when Ji — J3 > Jo -I- J2 and Ji — J3 < —(Jo + J2) 
(shaded regions in Fig. These midgap modes of Majo- 
rana fermions have two-fold degeneracy due to the time- 
reversal symmetry, or equivalently, spin 1/2 degrees of 
freedom bound on each edge. 



The Hamiltonian of one-dimensional system (105) has 
the same symmetry as that of two-dimensional system 
(p5|). Thus the I2 invariant is the product of 



sgn[$(q)] 



(108) 



Vi[w{q)] 

at the TRIM in the one-dimensional Brillouin zone, 

$(0) = Jo _ J, + J2 + J3, (109a) 



The phase boundaries obtained from the Z2 invariant 
coincides with those from numerics, which are alr eady 
given in ( |107| ). Since the class Dili Hamiltonian ( |105| ) 
can be decomposed into two independent blocks, each of 
which is a member of class AIII, the Z2 invariant in this 
case coincides with the even-odd parity of the integral 
invariant (winding number) of class AIII for the blocks. 
In turn, the winding number can be obtained by drawing 
the loop tr aject ory in the parameter space defined by 
$(q) inEq. (105) [see for discussion around Eq. ([72|)]; The 
product sgn[$(0)]sgn[$(7r)] (i.e., the Z2 invariant) then 
tells us that, when negative, the loop trajectory drawn 
by $(g) encloses the origin an odd number of times. 



VII. CONCLUSIONS 

In this paper, we have introduced a time-reversal sym- 
metric two-dimensional quantum spin model in a topo- 
logically non-trivial gapped phase, as a 7-matrix exten- 
sion of the Kitaev model on the square lattice. Through 
a fermion representation of the 7 matrices, this model is 
equivalent to a time-reversal symmetric two-dimensional 
topological superconductor (i.e., a system in class Dili 
in the Altland-Zirnbauer classification). The Hamilto- 
nian consists of nearest-neighbor interaction terms and 
next-nearest-neighbor interaction terms, all of which are 
transformed to free Majorana fermion hopping Hamil- 
tonian with Z2 gauge field. In the parameter space of 
the Hamiltonian, topologically trivial ground states and 
non-trivial ones are realized. We have confirmed that 
these two phases are distinguished by the Z2 topological 
invariant. 

We have shown using both numerical and analyti- 
cal methods the existence of topologically protected, a 
Kramers' pair of Majorana edge modes, which is a hall- 
mark of a time-reversal symmetric topological supercon- 
ductor. A local operator of a product of two flavors of 
Majorana operators, which is equivalent to the density 
operator of a complex fermion, has a nonvanishing corre- 
lation that decays in inverse-square of the distance along 
the edge and decays exponentially in the bulk. We have 
also shown numerically that a vortex of the Z2 gauge field 
hosts a Kramers' pair of zero-energy Majorana states. 

On the one-dimensional ladder lattice, we have con- 
structed the same type of extended Kitaev model. Sim- 
ilarly to the square lattice case, two topologically dis- 
tinct types of ground states appear in the phase diagram, 
which are characterized by the Z2 topological invariant. 
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Appendix: Jordan- Wigner transformation of the 
gamma matrix Kitaev model 

In this appendix, we present a solution to the Kitaev- 
type model in terms of the Jordan- Wigner trans- 
formation, following the solution of the QiLginal Kitaev 
model by Jordan- Wigner transformation. l3c3 



1. Jordan- Wigner transformation of the Dirac 
matrices 

As a first step, we note that it is possible to represent 
the Dirac matrices for a given site in terms of two complex 
fermions ci and C2 as 



(ci + 4), 

(ci - c\)/i, 



4), 



(C2 

(C2 - 4)/i, 



(A.l) 



where {ci, c\} = {c2, c\} = 1, and {ci, C2} =_[ci, 4} — 0. 
The right-hand side of four equations in (A.l) are Majo- 
rana fermion operators. Similarly, using C,^^ 



■"^(ci+cl), c' = 



S*'^"1(C2 - c\)/i, 

-e*™i(c2 + 4), 



(A.2) 



(A.3) 



where ni — 4ci, ?^2 — c\c2, and we have used the rela- 
tion e-"" = l-2na = -(ca + 4)(ca - 4). 

For the Kitaev- type model (0) , we need to prepare a set 
of two complex fermions cij and C2j for each site labeled 
by j. Accordingly, one needs to introduce string opera- 
tors to ensure commutation relations for operators sitting 
on different sites. At fist, we consider the string opera- 
tors for zeroth and second components of the Dirac ma- 
trices, which contain single Majorana operators made of 
ci. Since 0-links and 2-links horizontally connect neigh- 
boring sites (Fig. |l|), we define an order of sites on the 
two-dimensional lattice that runs horizontally: 



{xi,yi) < (x2,y2) <^ 

2/1 > y2 or (j/i = y2 and Xi < X2) 



(A.4) 



where x and y are the two-dimensional coordinates 
[Fig. |9|(a)]. Multiplying the right-hand side of the ze- 
roth and second components in Eqs. ( [A.l| ) and (|A.2|) by 
stringoperators of products of e*'^"^ with the order ( [A^ ) 
[Fig. 1(a)], 



t/, = n 

3>k 



(A.5) 



1'^ 



(a) 



(b) 



FIG. 9. String operators of Jordan- Wigner fermions that rep- 
resent Dirac matrices for (a) Majorana fermion made of ci 
(Uj), and (b) those made of C2 (Vj). 



makes these operators commute between different sites. 

Next, we can also introduce string operators for the 
first and third components, with another order of sites 
on the same two-dimensional lattice [Fig. ^(b)] 

(a;i,yi) < (2:2, 2/2) ^ 
xi < X2 or (xi — X2 and yi > 2/2), (^-6) 



as 



j>k 



(A.7) 



The latter string operator (A.7) ensures commutativity 
of the first and third components of Dirac matrix on dif- 
ferent sites. Finally, in order for Jordan- Wigner fermion 
representation of all four components of Dirac matrix to 
be bosonic, the first and third components of Dirac ma- 
trix on a site j is multiplied by string operator that runs 
all the sites besides j: 



(A., 



As a result, we obtain the Jordan- Wigner fermion repre- 
sentation of Dirac matrices as 



(Clj 
(C2j 



c. 

-c\-)Uj, 



-i{c2j 



(A.9) 



and 



(A.IO) 



cL-)e""'^V,-W^i 



Jordan- Wigner transformation of the 7 matrix 
Kitaev model 



By utilizing the Jordan- Wigner fermion representation 
of Dirac matrices, we can transform our spin model to a 
free Majorana fermion model, without redundancy. 
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On 0-links, since j € A and k € B co nnected by a 
0-link have the order j < k defined in (A.4), the nearest- 
neighbor interaction terms are transformed as follows. 



[Clk 



Consider ing t he or der o f j £ A and fc € B in the string 
operators ( (|A.4| ) and (|A.(]| )), nearest-neighbor interaction 
terms on the other three links are similarly transformed 
as 



cjcl 



2 2 



'{C2j - 4j)(c2fc -^4^), 



-(cij - 4j)(cifc + 4fc)' 



-(C2i - 4j)(c2fc +4fc)- 



Here we introduce two Majorana fermion operators on 
each site as 



(c2,+4^.)e*""^^ 

t 



-«(C2J - 4i)' 



C2fe 



(A.ll) 
(A.12) 



where j £ A and fc S B. Those Majorana fermion op- 
erators transform the nearest-neighbor interaction terms 



on 1-links and 3-links to free Majorana hopping terms 
without Z2 gauge operators, as 



/.i— 1,3 /j,-links 

= * E -^A. E (AlA^. + A^I). 

/X— 1,3 /i-links 



(A.13) 



Majorana fermions in Eq. ( A.12 ) also transform the 
nearest-neighbor interaction terms on 0-links and 2-links 
to free Majorana hopping terms, however in this case, 
with Z2 gauge operators: 

- E E + cj'cn 

/J,— 0,2 /j,-links 

= i E -^/^ E u,k{\',\i+\)\l). (A.14) 

/i— 0,2 /i-links 

where Ujk are Z2 gauge operators of the form 

Ujk = -(Cy + 4j)(c2j + 4j)(cifc ~ 4fc)(c2fc - 4fc)- 

(A.15) 

The Z2 gauge operators commute with each other, and 
\j, A|;, A^, A|: 



also commute with A|, At, A^, A|: 



[Ujk, Uirn] = 0, 

[ujkAt] = Kfc,Af] = 0. 



(A.16) 
(A.17) 



Consequently, we obtain the Jordan- Wigner transfor- 
mation of the 7 matrix Kitaev model as 



n 



E •^oWjfc+ E 

Vo-links l-linl( 



Jl 



X (A)Al. + A^^A^) 



E J2Ujk+ E -^3) 

ks 2-links 0-links / 

(A.I8) 



In this representation, the number of the Z2 gauge op- 
erators reduces to the half of that in the local Majorana 
representation ([l^), and therefore there is no redundancy 
in this fermionic representation. If we consider the peri- 
odic boundary condition, string operators appear in the 
Hamiltonian in connecting bath edges, as discussed in 
the case of the Kitaev model.l23 
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